#include <cstdlib>
#include <sys/types.h>
#include <cstdio>
#include <cmath>
#include <limits>

#include "polylib_templated.hh"
#define MAX_ITERATIONS  30
#define sign(a,b) ((b)<0 ? -fabs(a) : fabs(a))

#define _OPENMP
#ifdef _OPENMP
#define POLYLIB_USE_OPENMP
#endif

#ifdef POLYLIB_USE_OPENMP
#include <omp.h>
#define OPENMPMODE dynamic
#else
#define omp_get_num_threads() 0
#define omp_get_thread_num() 0
#endif


/*
LIBRARY ROUTINES FOR ORTHOGONAL POLYNOMIAL CALCULUS AND INTERPOLATION

  Spencer Sherwin
  Aeronautics, Imperial College London

  Based on codes by Einar Ronquist and Ron Henderson

  Abbreviations
  - z    -   Set of collocation/quadrature points
  - w    -   Set of quadrature weights
  - D    -   Derivative matrix
  - h    -   Lagrange Interpolant
  - I    -   Interpolation matrix
  - g    -   Gauss
  - gr   -   Gauss-Radau
  - gl   -   Gauss-Lobatto
  - j    -   Jacobi
  - l    -   Legendre  (Jacobi with alpha = beta =  0.0)
  - c    -   Chebychev (Jacobi with alpha = beta = -0.5)
  - m    -   point at minus 1 in Radau rules
  - p    -   point at plus  1 in Radau rules

  -----------------------------------------------------------------------
                        M A I N     R O U T I N E S
  -----------------------------------------------------------------------

  Points and Weights:

  zwgj        Compute Gauss-Jacobi         points and weights
  zwgrjm      Compute Gauss-Radau-Jacobi   points and weights (z=-1)
  zwgrjp      Compute Gauss-Radau-Jacobi   points and weights (z= 1)
  zwglj       Compute Gauss-Lobatto-Jacobi points and weights

  Derivative Matrices:

  Dgj         Compute Gauss-Jacobi         derivative matrix
  Dgrjm       Compute Gauss-Radau-Jacobi   derivative matrix (z=-1)
  Dgrjp       Compute Gauss-Radau-Jacobi   derivative matrix (z= 1)
  Dglj        Compute Gauss-Lobatto-Jacobi derivative matrix

  Lagrange Interpolants:

  hgj         Compute Gauss-Jacobi         Lagrange interpolants
  hgrjm       Compute Gauss-Radau-Jacobi   Lagrange interpolants (z=-1)
  hgrjp       Compute Gauss-Radau-Jacobi   Lagrange interpolants (z= 1)
  hglj        Compute Gauss-Lobatto-Jacobi Lagrange interpolants

  Interpolation Operators:

  Imgj        Compute interpolation operator gj->m
  Imgrjm      Compute interpolation operator grj->m (z=-1)
  Imgrjp      Compute interpolation operator grj->m (z= 1)
  Imglj       Compute interpolation operator glj->m

  Polynomial Evaluation:

  jacobfd     Returns value and derivative of Jacobi poly. at point z
  jacobd      Returns derivative of Jacobi poly. at point z (valid at z=-1,1)

  -----------------------------------------------------------------------
                    L O C A L      R O U T I N E S
  -----------------------------------------------------------------------

  jacobz      Returns Jacobi polynomial zeros
  gammaf      Gamma function for integer values and halves

  -----------------------------------------------------------------------
                        M A C R O S
  -----------------------------------------------------------------------

  Legendre  polynomial alpha = beta = 0
  Chebychev polynomial alpha = beta = -0.5

  Points and Weights:

  zwgl        Compute Gauss-Legendre          points and weights
  zwgrlm      Compute Gauss-Radau-Legendre    points and weights (z=-1)
  zwgrlp      Compute Gauss-Radau-Legendre    points and weights (z=+1)
  zwgll       Compute Gauss-Lobatto-Legendre  points and weights

  zwgc        Compute Gauss-Chebychev         points and weights
  zwgrcm      Compute Gauss-Radau-Chebychev   points and weights (z=-1)
  zwgrcp      Compute Gauss-Radau-Chebychev   points and weights (z=+1)
  zwglc       Compute Gauss-Lobatto-Chebychev points and weights

  Derivative Operators:

  Dgl         Compute Gauss-Legendre          derivative matrix
  Dgrlm       Compute Gauss-Radau-Legendre    derivative matrix (z=-1)
  Dgrlp       Compute Gauss-Radau-Legendre    derivative matrix (z=+1)
  Dgll        Compute Gauss-Lobatto-Legendre  derivative matrix

  Dgc         Compute Gauss-Chebychev         derivative matrix
  Dgrcm       Compute Gauss-Radau-Chebychev   derivative matrix (z=-1)
  Dgrcp       Compute Gauss-Radau-Chebychev   derivative matrix (z=+1)
  Dglc        Compute Gauss-Lobatto-Chebychev derivative matrix

  Lagrangian Interpolants:

  hgl         Compute Gauss-Legendre          Lagrange interpolants
  hgrlm       Compute Gauss-Radau-Legendre    Lagrange interpolants (z=-1)
  hgrlp       Compute Gauss-Radau-Legendre    Lagrange interpolants (z=+1)
  hgll        Compute Gauss-Lobatto-Legendre  Lagrange interpolants

  hgc         Compute Gauss-Chebychev         Lagrange interpolants
  hgrcm       Compute Gauss-Radau-Chebychev   Lagrange interpolants (z=-1)
  hgrcp       Compute Gauss-Radau-Chebychev   Lagrange interpolants (z=+1)
  hglc        Compute Gauss-Lobatto-Chebychev Lagrange interpolants

  Interpolation Operators:

  Imgl        Compute interpolation operator gl->m
  Imgrlm      Compute interpolation operator grl->m (z=-1)
  Imgrlp      Compute interpolation operator grl->m (z=+1)
  Imgll       Compute interpolation operator gll->m

  Imgc        Compute interpolation operator gc->m
  Imgrcm      Compute interpolation operator grc->m (z=-1)
  Imgrcp      Compute interpolation operator grc->m (z=+1)
  Imglc       Compute interpolation operator glc->m

  ------------------------------------------------------------------------

  Useful references:

  - [1] Gabor Szego: Orthogonal Polynomials, American Mathematical Society,
      Providence, Rhode Island, 1939.
  - [2] Abramowitz \& Stegun: Handbook of Mathematical Functions,
      Dover, New York, 1972.
  - [3] Canuto, Hussaini, Quarteroni \& Zang: Spectral Methods in Fluid
      Dynamics, Springer-Verlag, 1988.
  - [4] Ghizzetti \& Ossicini: Quadrature Formulae, Academic Press, 1970.
  - [5] Karniadakis \& Sherwin: Spectral/hp element methods for CFD, 1999


  NOTES
  -----
  (1) All routines are double precision.
  (2) All array subscripts start from zero, i.e. vector[0..N-1]
*/


namespace polylib {


  /**
  \brief  Gauss-Jacobi zeros and weights.

  \li Generate \a np Gauss Jacobi zeros, \a z, and weights,\a w,
  associated with the Jacobi polynomial \f$ P^{\alpha,\beta}_{np}(z)
  \f$,

  \li Exact for polynomials of order \a 2np-1 or less
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::zwgj (Jacobi<BaseType>::UnmanagedArray z, 
                                        Jacobi<BaseType>::UnmanagedArray w,
                                        const unsigned int np, 
                                        const BaseType alpha, const BaseType beta)
  {
    register int i;
    BaseType   fac;
    const BaseType   one = 1.0, apb = alpha + beta, two = 2.0;

    jacobz (np, z, alpha, beta);
    jacobd (np, z, w, np, alpha, beta);

    fac  = pow(two,apb + one)*Utility::gammaF(alpha + np + one)*Utility::gammaF(beta + np + one);
    fac /= Utility::gammaF(np + one)*Utility::gammaF(apb + np + one);

    for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]*(one-z[i]*z[i]));

    return;
  }


  /**
  \brief  Gauss-Radau-Jacobi zeros and weights with end point at \a z=-1.

  \li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
  associated with the  polynomial \f$(1+z) P^{\alpha,\beta+1}_{np-1}(z)
  \f$.

  \li  Exact for polynomials of order \a 2np-2 or less
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::zwgrjm(Jacobi<BaseType>::UnmanagedArray z, 
                                        Jacobi<BaseType>::UnmanagedArray w,
                                        const unsigned int np, 
                                        const BaseType alpha, const BaseType beta)
  {
    if(np == 1){
      z[0] = 0.0;
      w[0] = 2.0;
    }
    else{
      register int i;
      BaseType   fac;
      const BaseType   one = 1.0, apb = alpha + beta, two = 2.0;

      z[0] = -one;
      jacobz  (np-1, &z[1], alpha, beta+1.);
      jacobfd (np, z, &w[0], NULL, np-1, alpha, beta);

      fac  = pow(two,apb)*Utility::gammaF(alpha + np)*Utility::gammaF(beta + np);
      fac /= Utility::gammaF(np*1.)*(beta + np)*Utility::gammaF(apb + np + 1);

      for(i = 0; i < np; ++i) w[i] = fac*(1-z[i])/(w[i]*w[i]);
      w[0] *= (beta + one);
    }

    return;
  }


  /**
  \brief  Gauss-Radau-Jacobi zeros and weights with end point at \a z=1


  \li Generate \a np Gauss-Radau-Jacobi zeros, \a z, and weights,\a w,
  associated with the  polynomial \f$(1-z) P^{\alpha+1,\beta}_{np-1}(z)
  \f$.

  \li Exact for polynomials of order \a 2np-2 or less
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::zwgrjp(Jacobi<BaseType>::UnmanagedArray z, 
                                        Jacobi<BaseType>::UnmanagedArray w,
                                        const unsigned int np, 
                                        const BaseType alpha, const BaseType beta)
  {
    if(np == 1){
      z[0] = 0.0;
      w[0] = 2.0;
    }
    else{
      register int i;
      BaseType   fac;
      const BaseType   one = 1.0, apb = alpha + beta, two = 2.0;

      jacobz  (np-1, &z[0], alpha+1., beta);
      z[np-1] = one;
      jacobfd (np, z, &w[0], NULL, np-1, alpha, beta);

      fac  = pow(two,apb)*Utility::gammaF(alpha + np)*Utility::gammaF(beta + np);
      fac /= Utility::gammaF(np*1.)*(alpha + np)*Utility::gammaF(apb + np + 1);

      for(i = 0; i < np; ++i) w[i] = fac*(1+z[i])/(w[i]*w[i]);
      w[np-1] *= (alpha + one);
    }

    return;
  }


  /**
  \brief  Gauss-Lobatto-Jacobi zeros and weights with end point at \a z=-1,\a 1


  \li Generate \a np Gauss-Lobatto-Jacobi points, \a z, and weights, \a w,
  associated with polynomial \f$ (1-z)(1+z) P^{\alpha+1,\beta+1}_{np-2}(z) \f$
  \li Exact for polynomials of order \a 2np-3 or less
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::zwglj(Jacobi<BaseType>::UnmanagedArray z, 
                                        Jacobi<BaseType>::UnmanagedArray w,
                                        const unsigned int np, 
                                        const BaseType alpha, const BaseType beta)
  {
    if( np == 1 ){
      z[0] = 0.0;
      w[0] = 2.0;
    }
    else{
      register int i;
      BaseType   fac;
      const BaseType   one = 1.0, apb = alpha + beta, two = 2.0;

      z[0]    = -one;
      z[np-1] =  one;
      jacobz  (np-2, &z[1], alpha + one, beta + one);
      jacobfd (np, z, w, NULL, np-1, alpha, beta);

      fac  = pow(two,apb + 1)*Utility::gammaF(alpha + np)*Utility::gammaF(beta + np);
      fac /= (np-1)*Utility::gammaF(np*1.)*Utility::gammaF(alpha + beta + np + one);

      for(i = 0; i < np; ++i) w[i] = fac/(w[i]*w[i]);
      w[0]    *= (beta  + one);
      w[np-1] *= (alpha + one);
    }

    return;
  }


  /**
    \brief Compute the Derivative Matrix and its transpose associated
    with the Gauss-Jacobi zeros.

    \li Compute the derivative matrix, \a d, and its transpose, \a dt,
    associated with the n_th order Lagrangian interpolants through the
    \a np Gauss-Jacobi points \a z such that \n
    \f$  \frac{du}{dz}(z[i]) =  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$

    \li d and dt are both square matrices.
*/
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::Dgj(Jacobi<BaseType>::UnmanagedArray D, 
                                      Jacobi<BaseType>::UnmanagedArray Dt, 
                                      Jacobi<BaseType>::UnmanagedArray z, 
                                      const unsigned int np, 
                                      const BaseType m_alpha, const BaseType m_beta)
  {
    const BaseType one = 1.0, two = 2.0;

    if (np <= 0){
      D[0] = Dt[0] = 0.0;
    }
    else{
      register int i,j;

      Jacobi<BaseType>::UnmanagedArray pd;
      pd = (Jacobi<BaseType>::UnmanagedArray)malloc(np*sizeof(BaseType)) ;

      jacobd(np,z,pd,np,m_alpha,m_beta);

      for (i = 0; i < np; i++){
        for (j = 0; j < np; j++){

          if (i != j)
          D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
          else
          D[i*np+j] = (m_alpha - m_beta + (m_alpha + m_beta + two)*z[i])/
          (two*(one - z[i]*z[i]));

          Dt[j*np+i] = D[i*np+j];
        }
      }
      free(pd);
    }
    return;
  }


  /**
    \brief Compute the Derivative Matrix and its transpose associated
    with the Gauss-Radau-Jacobi zeros with a zero at \a z=-1.

    \li Compute the derivative matrix, \a d, and its transpose, \a dt,
    associated with the n_th order Lagrangian interpolants through the
    \a np Gauss-Radau-Jacobi points \a z such that \n
    \f$  \frac{du}{dz}(z[i]) =  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$

    \li d and dt are both square matrices.
*/
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::Dgrjm(Jacobi<BaseType>::UnmanagedArray D, 
                                        Jacobi<BaseType>::UnmanagedArray Dt, 
                                        Jacobi<BaseType>::UnmanagedArray z, 
                                        const unsigned int np, 
                                        const BaseType m_alpha, const BaseType m_beta)
  {
    if (np <= 0){
      D[0] = Dt[0] = 0.0;
    }
    else{
      register int i, j;
      const BaseType   one = 1.0, two = 2.0;
      Jacobi<BaseType>::UnmanagedArray pd;
      pd = (Jacobi<BaseType>::UnmanagedArray)malloc(np*sizeof(BaseType)) ;

      pd[0] = pow(-one,np-1)*Utility::gammaF(np+m_beta+one);
      pd[0] /= Utility::gammaF(np*1.)*Utility::gammaF(m_beta+two);
      jacobd(np-1,z+1,pd+1,np-1,m_alpha,m_beta+1);
      for(i = 1; i < np; ++i) pd[i] *= (1+z[i]);

      for (i = 0; i < np; i++)
      for (j = 0; j < np; j++){
        if (i != j)
        D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
        else {
          if(i == 0)
          D[i*np+j] = -(np + m_alpha + m_beta + one)*(np - one)/
          (two*(m_beta + two));
          else
          D[i*np+j] = (m_alpha - m_beta + one + (m_alpha + m_beta + one)*z[i])/
          (two*(one - z[i]*z[i]));
        }

        Dt[j*np+i] = D[i*np+j];
      }
      free(pd);
    }

    return;
  }


  /**
    \brief Compute the Derivative Matrix and its transpose associated
    with the Gauss-Radau-Jacobi zeros with a zero at \a z=1.

    \li Compute the derivative matrix, \a d, and its transpose, \a dt,
    associated with the n_th order Lagrangian interpolants through the
    \a np Gauss-Radau-Jacobi points \a z such that \n
    \f$  \frac{du}{dz}(z[i]) =  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$

    \li d and dt are both square matrices.
*/
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::Dgrjp(Jacobi<BaseType>::UnmanagedArray D, 
                                        Jacobi<BaseType>::UnmanagedArray Dt, 
                                        Jacobi<BaseType>::UnmanagedArray z, 
                                        const unsigned int np, 
                                        const BaseType m_alpha, const BaseType m_beta)
  {
    if (np <= 0){
      D[0] = Dt[0] = 0.0;
    }
    else{
      register int i, j;
      const BaseType   one = 1.0, two = 2.0;
      Jacobi<BaseType>::UnmanagedArray pd;
      pd = (Jacobi<BaseType>::UnmanagedArray)malloc(np*sizeof(BaseType)) ;

      jacobd(np-1,z,pd,np-1,m_alpha+1,m_beta);
      for(i = 0; i < np-1; ++i) pd[i] *= (1-z[i]);
      pd[np-1] = -Utility::gammaF(np+m_alpha+one);
      pd[np-1] /= Utility::gammaF(np*1.)*Utility::gammaF(m_alpha+two);

      for (i = 0; i < np; i++)
      for (j = 0; j < np; j++){
        if (i != j)
        D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
        else {
          if(i == np-1)
          D[i*np+j] = (np + m_alpha + m_beta + one)*(np - one)/
          (two*(m_alpha + two));
          else
          D[i*np+j] = (m_alpha - m_beta - one + (m_alpha + m_beta + one)*z[i])/
          (two*(one - z[i]*z[i]));
        }

        Dt[j*np+i] = D[i*np+j];
      }
      free(pd);
    }

    return;
  }

  /**
    \brief Compute the Derivative Matrix and its transpose associated
    with the Gauss-Lobatto-Jacobi zeros.

    \li Compute the derivative matrix, \a d, and its transpose, \a dt,
    associated with the n_th order Lagrangian interpolants through the
    \a np Gauss-Lobatto-Jacobi points \a z such that \n
    \f$  \frac{du}{dz}(z[i]) =  \sum_{j=0}^{np-1} D[i*np+j] u(z[j]) \f$

    \li d and dt are both square matrices.
*/
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::Dglj( Jacobi<BaseType>::UnmanagedArray  D, 
                                        Jacobi<BaseType>::UnmanagedArray  Dt, 
                                        Jacobi<BaseType>::UnmanagedArray z, 
                                        const unsigned int np, 
                                        const BaseType m_alpha, const BaseType m_beta )
  {
    if (np <= 0){
      D[0] = Dt[0] = 0.0;
    }
    else{
      register int i, j;
      const BaseType   one = 1.0, two = 2.0;
      Jacobi<BaseType>::UnmanagedArray pd;
      pd = (Jacobi<BaseType>::UnmanagedArray)malloc(np*sizeof(BaseType)) ;

      pd[0]  = two*pow(-one,np)*Utility::gammaF(np + m_beta);
      pd[0] /= Utility::gammaF(np - one)*Utility::gammaF(m_beta + two);
      jacobd(np-2,z+1,pd+1,np-2,m_alpha+1,m_beta+1);
      for(i = 1; i < np-1; ++i) pd[i] *= (one-z[i]*z[i]);
      pd[np-1]  = -two*Utility::gammaF(np + m_alpha);
      pd[np-1] /= Utility::gammaF(np - one)*Utility::gammaF(m_alpha + two);

      for (i = 0; i < np; i++)
      for (j = 0; j < np; j++){
        if (i != j)
        D[i*np+j] = pd[i]/(pd[j]*(z[i]-z[j]));
        else {
          if      (i == 0)
          D[i*np+j] = (m_alpha - (np-1)*(np + m_alpha + m_beta))/(two*(m_beta+ two));
          else if (i == np-1)
          D[i*np+j] =-(m_beta - (np-1)*(np + m_alpha + m_beta))/(two*(m_alpha+ two));
          else
          D[i*np+j] = (m_alpha - m_beta + (m_alpha + m_beta)*z[i])/
          (two*(one - z[i]*z[i]));
        }

        Dt[j*np+i] = D[i*np+j];
      }
      free(pd);
    }

    return;
  }


  /**
    \brief Compute the value of the \a i th Lagrangian interpolant through
    the \a np Gauss-Jacobi points \a zgj at the arbitrary location \a z.

    \li \f$ -1 \leq z \leq 1 \f$

    \li Uses the defintion of the Lagrangian interpolant:\n
    %
    \f$ \begin{array}{rcl}
    h_j(z) =  \left\{ \begin{array}{ll}
    \displaystyle \frac{P_{np}^{\m_alpha,\beta}(z)}
    {[P_{np}^{\m_alpha,\beta}(z_j)]^\prime
    (z-z_j)} & \mbox{if $z \ne z_j$}\\
    & \\
    1 & \mbox{if $z=z_j$}
    \end{array}
    \right.
    \end{array}   \f$
  */
  template<typename BaseType>
  BaseType Jacobi<BaseType>::Utility::hgj ( BaseType z, 
                                            BaseType  zi, 
                                            const unsigned int np, 
                                            const BaseType m_alpha, 
                                            const BaseType m_beta)
  {
    BaseType dz, p, pd, h;

    dz  = z - zi;
    if (fabs(dz) < std::numeric_limits<BaseType>::epsilon()*100) return 1.0;

    jacobd (1, &zi, &pd , np, m_alpha, m_beta);
    jacobfd(1, &z , &p, NULL , np, m_alpha, m_beta);
    h = p/(pd*dz);

    return h;
  }

  /**
    \brief Compute the value of the \a i th Lagrangian interpolant through the
    \a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
    \a z. This routine assumes \a zgrj includes the point \a -1.

    \li \f$ -1 \leq z \leq 1 \f$

    \li Uses the defintion of the Lagrangian interpolant:\n
    %
    \f$ \begin{array}{rcl}
    h_j(z) = \left\{ \begin{array}{ll}
    \displaystyle \frac{(1+z) P_{np-1}^{\alpha,\beta+1}(z)}
    {((1+z_j) [P_{np-1}^{\alpha,\beta+1}(z_j)]^\prime +
    P_{np-1}^{\alpha,\beta+1}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
    & \\
    1 & \mbox{if $z=z_j$}
    \end{array}
    \right.
    \end{array}   \f$
  */
  template<typename BaseType>
  BaseType Jacobi<BaseType>::Utility::hgrjm ( BaseType z, 
                                              BaseType zi, 
                                              const unsigned int np, 
                                              const BaseType alpha, 
                                              const BaseType beta)
  {
    BaseType p, pd, h;

    BaseType dz  = z - zi;
    if (fabs(dz) < std::numeric_limits<BaseType>::epsilon()*100) return 1.0;

    jacobfd (1, &zi, &p , NULL, np-1, alpha, beta + 1);
    // need to use this routine in case zi = -1 or 1
    jacobd  (1, &zi, &pd, np-1, alpha, beta + 1);
    h = (1.0 + zi)*pd + p;
    jacobfd (1, &z, &p, NULL,  np-1, alpha, beta + 1);
    h = (1.0 + z )*p/(h*dz);

    return h;
  }


  /**
    \brief Compute the value of the \a i th Lagrangian interpolant through the
    \a np Gauss-Radau-Jacobi points \a zgrj at the arbitrary location
    \a z. This routine assumes \a zgrj includes the point \a +1.

    \li \f$ -1 \leq z \leq 1 \f$

    \li Uses the defintion of the Lagrangian interpolant:\n
%
    \f$ \begin{array}{rcl}
    h_j(z) = \left\{ \begin{array}{ll}
    \displaystyle \frac{(1-z) P_{np-1}^{\alpha+1,\beta}(z)}
    {((1-z_j) [P_{np-1}^{\alpha+1,\beta}(z_j)]^\prime -
    P_{np-1}^{\alpha+1,\beta}(z_j) ) (z-z_j)} & \mbox{if $z \ne z_j$}\\
    & \\
    1 & \mbox{if $z=z_j$}
    \end{array}
    \right.
    \end{array}   \f$
*/

  template<typename BaseType>
  BaseType Jacobi<BaseType>::Utility::hgrjp (BaseType z, 
                                             BaseType  zi, 
                                             const unsigned int np, 
                                             const BaseType m_alpha, 
                                             const BaseType m_beta)
  {
    BaseType dz, p, pd, h;

    //zi  = zgrj[i];
    dz  = z - zi;
    if (fabs(dz) < std::numeric_limits<BaseType>::epsilon()*100) return 1.0;

    jacobfd (1, &zi, &p , NULL, np-1, m_alpha+1, m_beta );
    // need to use this routine in case z = -1 or 1
    jacobd  (1, &zi, &pd, np-1, m_alpha+1, m_beta );
    h = (1.0 - zi)*pd - p;
    jacobfd (1, &z, &p, NULL,  np-1, m_alpha+1, m_beta);
    h = (1.0 - z )*p/(h*dz);

    return h;
  }


  /**
    \brief Compute the value of the \a i th Lagrangian interpolant through the
    \a np Gauss-Lobatto-Jacobi points \a zgrj at the arbitrary location
    \a z.

    \li \f$ -1 \leq z \leq 1 \f$

    \li Uses the defintion of the Lagrangian interpolant:\n
%
    \f$ \begin{array}{rcl}
    h_j(z) = \left\{ \begin{array}{ll}
    \displaystyle \frac{(1-z^2) P_{np-2}^{\alpha+1,\beta+1}(z)}
    {((1-z^2_j) [P_{np-2}^{\alpha+1,\beta+1}(z_j)]^\prime -
    2 z_j P_{np-2}^{\alpha+1,\beta+1}(z_j) ) (z-z_j)}&\mbox{if $z \ne z_j$}\\
    & \\
    1 & \mbox{if $z=z_j$}
    \end{array}
    \right.
    \end{array}   \f$
*/

  template<typename BaseType>
  BaseType Jacobi<BaseType>::Utility::hglj (BaseType z, 
                                            BaseType zi, 
                                            const unsigned int np, 
                                            const BaseType m_alpha, 
                                            const BaseType m_beta)
  {
    const BaseType one = 1., two = 2.;
    BaseType dz, p, pd, h;

    //  zi  = zglj[i];
    dz  = z - zi;
    if (fabs(dz) < std::numeric_limits<BaseType>::epsilon()*100) return 1.0;

    jacobfd(1, &zi, &p , NULL, np-2, m_alpha + one, m_beta + one);
    // need to use this routine in case z = -1 or 1
    jacobd (1, &zi, &pd, np-2, m_alpha + one, m_beta + one);
    h = (one - zi*zi)*pd - two*zi*p;
    jacobfd(1, &z, &p, NULL, np-2, m_alpha + one, m_beta + one);
    h = (one - z*z)*p/(h*dz);

    return h;
  }


  /**
    \brief Routine to calculate Jacobi polynomials, \f$
    P^{\alpha,\beta}_n(z) \f$, and their first derivative, \f$
    \frac{d}{dz} P^{\alpha,\beta}_n(z) \f$.

    \li This function returns the vectors \a poly_in and \a poly_d
    containing the value of the \f$ n^th \f$ order Jacobi polynomial
    \f$ P^{\alpha,\beta}_n(z) \alpha > -1, \beta > -1 \f$ and its
    derivative at the \a np points in \a z[i]

    - If \a poly_in = NULL then only calculate derivatice

    - If \a polyd   = NULL then only calculate polynomial

    - To calculate the polynomial this routine uses the recursion
    relationship (see appendix A ref [4]) :
    \f$ \begin{array}{rcl}
    P^{\alpha,\beta}_0(z) &=& 1 \\
    P^{\alpha,\beta}_1(z) &=& \frac{1}{2} [ \alpha-\beta+(\alpha+\beta+2)z] \\
    a^1_n P^{\alpha,\beta}_{n+1}(z) &=& (a^2_n + a^3_n z)
    P^{\alpha,\beta}_n(z) - a^4_n P^{\alpha,\beta}_{n-1}(z) \\
    a^1_n &=& 2(n+1)(n+\alpha + \beta + 1)(2n + \alpha + \beta) \\
    a^2_n &=& (2n + \alpha + \beta + 1)(\alpha^2 - \beta^2)  \\
    a^3_n &=& (2n + \alpha + \beta)(2n + \alpha + \beta + 1)
    (2n + \alpha + \beta + 2)  \\
    a^4_n &=& 2(n+\alpha)(n+\beta)(2n + \alpha + \beta + 2)
    \end{array} \f$

    - To calculate the derivative of the polynomial this routine uses
    the relationship (see appendix A ref [4]) :
    \f$ \begin{array}{rcl}
    b^1_n(z)\frac{d}{dz} P^{\alpha,\beta}_n(z)&=&b^2_n(z)P^{\alpha,\beta}_n(z)
    + b^3_n(z) P^{\alpha,\beta}_{n-1}(z) \hspace{2.2cm} \\
    b^1_n(z) &=& (2n+\alpha + \beta)(1-z^2) \\
    b^2_n(z) &=& n[\alpha - \beta - (2n+\alpha + \beta)z]\\
    b^3_n(z) &=& 2(n+\alpha)(n+\beta)
    \end{array} \f$

    - Note the derivative from this routine is only valid for -1 < \a z < 1.
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::JacobFunctionDerivative(  const int np, 
                                                            const Jacobi<BaseType>::UnmanagedArray  z, 
                                                            Jacobi<BaseType>::UnmanagedArray poly_in, 
                                                            Jacobi<BaseType>::UnmanagedArray polyd, 
                                                            const unsigned int n, 
                                                            const BaseType alpha, 
                                                            const BaseType beta)
  {
    const BaseType  zero = 0.0, one = 1.0, two = 2.0, apb = alpha + beta;

    if(!np)
      return;

    if(n == 0){
      if(poly_in)
      {
#ifdef POLYLIB_USE_OPENMP
#pragma omp parallel for shared(poly_in) schedule(OPENMPMODE, 1) // if (np>2)
#endif
        for(int i = 0; i < np; ++i)
          poly_in[i] = 1.0;
      }
      if(polyd)
      {
#ifdef POLYLIB_USE_OPENMP
#pragma omp parallel for shared(polyd) schedule(OPENMPMODE, 1) // if (np>2)
#endif
        for(int i = 0; i < np; ++i)
          polyd[i] = 0.0;
      }
    }
    else if (n == 1){
      if(poly_in)
      {
#ifdef POLYLIB_USE_OPENMP
#pragma omp parallel for shared(poly_in, z) schedule(OPENMPMODE, 1) //  if (np>2)
#endif
        for(int i = 0; i < np; ++i)
          poly_in[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
      }
      if(polyd)
      {
#ifdef POLYLIB_USE_OPENMP
#pragma omp parallel for shared(polyd) schedule(OPENMPMODE, 1) // if (np>2)
#endif
        for(int i = 0; i < np; ++i)
          polyd[i] = 0.5*(alpha + beta + two);
      }
    }
    else{

      BaseType   *poly, *polyn1; //,*polyn2;

      if(poly_in)
      {
        // switch for case of no poynomial function return
        polyn1 = (Jacobi<BaseType>::UnmanagedArray)malloc(2*np*sizeof(double));
        //polyn2 = polyn1+np;
        poly   = poly_in;
      }
      else
      {
        polyn1 = (Jacobi<BaseType>::UnmanagedArray)malloc(3*np*sizeof(double));
        //polyn2 = polyn1+np;
        poly   = polyn1+2*np;
      }

#ifdef POLYLIB_USE_OPENMP
      //#pragma omp parallel shared(n, z, poly, polyn1, polyn2, polyd) if (np>2)
#endif
      {
#ifdef POLYLIB_USE_OPENMP
#pragma omp parallel for schedule(OPENMPMODE,1) shared(polyn1) firstprivate(z) // if (np>4)
#endif
        for(int i = 0; i < np; ++i) {
#ifdef POLYLIB_USE_OPENMP
#pragma omp critical
#endif
          {
            polyn1[i] = 0.5*(alpha - beta + (alpha + beta + two)*z[i]);
            polyn1[i+np] = one;
          }
        }

#ifdef POLYLIB_USE_OPENMP
        //#pragma omp flush(polyn1)
        //#pragma omp barrier
        //#pragma omp parallel for schedule(static, 2) shared(poly, np) firstprivate(z, polyn1) if (np>4)
#endif
        for(int k = 2; k <= n; ++k) {

          const BaseType a1 =  two*k*(k + apb)*(two*k + apb - two);
          const BaseType a2 = (two*k + apb - one)*(alpha*alpha - beta*beta) / a1 ;
          const BaseType a3 = (two*k + apb - two)*(two*k + apb - one)*(two*k + apb) / a1 ;
          const BaseType a4 =  two*(k + alpha - one)*(k + beta - one)*(two*k + apb) / a1 ;
          
          Jacobi<BaseType>::UnmanagedArray polyn2 = polyn1+np;

#ifdef POLYLIB_USE_OPENMP
          //#pragma omp for // parallel shared(np, z, polyn1, poly) private(a1, a2, a3, a4)
#pragma omp parallel for schedule(OPENMPMODE,1) shared(poly) firstprivate(z, polyn1, polyn2) // if (np>4)
#endif
          for(int i = 0; i < np; ++i) {

            BaseType temp = ( (a2 + a3*z[i])*polyn1[i] - a4*polyn2[i]) ;

#ifdef POLYLIB_USE_OPENMP
            //#pragma omp critical
#endif
            {
              poly[i] = temp;
              polyn2[i] = polyn1[i];
              polyn1[i] = temp;
              
            }
          }
          // flush all the data to memory
#ifdef POLYLIB_USE_OPENMP
#pragma omp flush(poly, polyn1, polyn2)
#endif
        }

        if(polyd) {

          const BaseType a1 = n*(alpha - beta) ;
          const BaseType a2 = n*(two*n + alpha + beta) / a1 ;
          const BaseType a3 = two*(n + alpha)*(n + beta) / a1 ;
          const BaseType a4 = (two*n + alpha + beta) / a1 ;

          // note polyn2 points to polyn1 at end of poly iterations
#ifdef POLYLIB_USE_OPENMP
#pragma omp parallel for shared(polyd) firstprivate(z, poly, polyn1) // if (np>4)
#endif
          for(int i = 0; i < np; ++i)
          {
            BaseType temp = (a1 - a2*z[i])*poly[i] + a3*polyn1[i+np] ;
            
#ifdef POLYLIB_USE_OPENMP
            #pragma omp critical
#endif
            polyd[i] = temp / (one - z[i]*z[i]) ;
          }
#ifdef POLYLIB_USE_OPENMP
#pragma omp flush(polyd)
#endif
        }

      } // end OMP pragma
      
      free(polyn1);
    } // else block for n==0

    return;
  }


  template<typename BaseType>
  void Jacobi<BaseType>::Utility::JacobFunction(const unsigned int np, 
                                                const Jacobi<BaseType>::UnmanagedArray  z, 
                                                Jacobi<BaseType>::UnmanagedArray  poly, 
                                                const unsigned int n, 
                                                const BaseType alpha, 
                                                const BaseType beta)
  {
    jacobfd(np,z,poly,NULL,n,alpha,beta);
  }

  /**
  \brief Calculate the  derivative of Jacobi polynomials

  \li Generates a vector \a poly of values of the derivative of the
  \a n th order Jacobi polynomial \f$ P^(\alpha,\beta)_n(z)\f$ at the
  \a np points \a z.

  \li To do this we have used the relation
  \n
  \f$ \frac{d}{dz} P^{\alpha,\beta}_n(z)
  = \frac{1}{2} (\alpha + \beta + n + 1)  P^{\alpha,\beta}_n(z) \f$

  \li This formulation is valid for \f$ -1 \leq z \leq 1 \f$
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::JacobDerivative(const unsigned int np, 
                                                  const Jacobi<BaseType>::UnmanagedArray z, 
                                                  Jacobi<BaseType>::UnmanagedArray polyd, 
                                                  const unsigned int n, 
                                                  const BaseType alpha, 
                                                  const BaseType beta)
  {
    register int i;
    const BaseType one = 1.0;
    if(n == 0)
    for(i = 0; i < np; ++i) polyd[i] = 0.0;
    else{
      jacobf(np,z,polyd,n-1,alpha+one,beta+one);
      //    jacobfd(np,z,polyd,NULL,n-1,alpha+one,beta+one);
      for(i = 0; i < np; ++i)
      polyd[i] *= 0.5*(alpha + beta + (double)n + one);
    }
    return;
  }


  /**
  \brief Calculate the Gamma function , \f$ \Gamma(n)\f$, for integer
  values and halves.

  Determine the value of \f$\Gamma(n)\f$ using:

  \f$ \Gamma(n) = (n-1)!  \mbox{ or  }  \Gamma(n+1/2) = (n-1/2)\Gamma(n-1/2)\f$

  where \f$ \Gamma(1/2) = \sqrt(\pi)\f$
  */
  template<typename BaseType>
  BaseType Jacobi<BaseType>::Utility::gammaF(BaseType x)
  {
    BaseType gamma = 1.0;

    if     (x == -0.5) gamma = -2.0*sqrt(M_PI);
    else if (!x) return gamma;
    else if ((x-(int)x) == 0.5){
      int n = (int) x;
      BaseType tmp = x;

      gamma = sqrt(M_PI);
      while(n--){
        tmp   -= 1.0;
        gamma *= tmp;
      }
    }
    else if ((x-(int)x) == 0.0){
      int n = (int) x;
      double tmp = x;

      while(--n){
        tmp   -= 1.0;
        gamma *= tmp;
      }
    }
    else
    fprintf(stderr,"%f is not of integer or half order\n",x);
    return gamma;
  }

  /**
    \brief  Calculate the \a n zeros, \a z, of the Jacobi polynomial, i.e.
    \f$ P_n^{\alpha,\beta}(z) = 0 \f$

    This routine is only value for \f$( \alpha > -1, \beta > -1)\f$
    and uses polynomial deflation in a Newton iteration
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::JacZerosNewton( const unsigned int n, 
                                                  Jacobi<BaseType>::UnmanagedArray z, 
                                                  BaseType alpha, BaseType beta)
  {
    if(n == 0)
      return;

    const BaseType   dth = M_PI/(2.0*(double)n);
    const BaseType one = 1.0, two = 2.0;

    BaseType   poly,pder,rlast=0.0;
    BaseType   delr,r;
    
    for(unsigned int k = 0; k < n; ++k){
      r = -cos((two*(double)k + one) * dth);
      if(k) r = 0.5*(r + rlast);

      for(unsigned int j = 1; j < MAX_ITERATIONS; ++j){
        jacobfd(1,&r,&poly, &pder, n, alpha, beta);

        BaseType sum = 0.0 ;
        for(unsigned i = 0; i < k; ++i)
        sum += one/(r - z[i]);

        delr = -poly / (pder - sum * poly);
        r   += delr;
        if( fabs(delr) < std::numeric_limits<BaseType>::epsilon()*100 ) break;
      }
      z[k]  = r;
      rlast = r;
    }
    return;
  }


  /**
  \brief Zero determination through the eigenvalues of a tridiagonal
  matrix from teh three term recursion relationship.

  Set up a symmetric tridiagonal matrix

  \f$ \left [  \begin{array}{ccccc}
  a[0] & b[0]   &        &        & \\
  b[0] & a[1]   & b[1]   &        & \\
    0   & \ddots & \ddots & \ddots &  \\
        &        & \ddots & \ddots & b[n-2] \\
        &        &        & b[n-2] & a[n-1] \end{array} \right ] \f$

  Where the coefficients a[n], b[n] come from the  recurrence relation

  \f$  b_j p_j(z) = (z - a_j ) p_{j-1}(z) - b_{j-1}   p_{j-2}(z) \f$

  where \f$ j=n+1\f$ and \f$p_j(z)\f$ are the Jacobi (normalized)
  orthogonal polynomials \f$ \alpha,\beta > -1\f$( integer values and
  halves). Since the polynomials are orthonormalized, the tridiagonal
  matrix is guaranteed to be symmetric. The eigenvalues of this
  matrix are the zeros of the Jacobi polynomial.
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::JacZerosQL(const unsigned int n, Jacobi<BaseType>::UnmanagedArray a, BaseType alpha, BaseType beta){
    
    if(n == 0)
      return;

    Jacobi<BaseType>::ManagedArray b (n) ;

    // generate normalised terms
    const BaseType apb  = alpha + beta;
    BaseType apbi = 2.0 + apb;
    const BaseType a2b2 = beta*beta-alpha*alpha;

    b[n-1] = pow(2.0,apb+1.0)*gammaF(alpha+1.0)*gammaF(beta+1.0)/gammaF(apbi);
    a[0]   = (beta-alpha)/apbi;
    b[0]   = sqrt(4.0*(1.0+alpha)*(1.0+beta)/((apbi+1.0)*apbi*apbi));

    for(unsigned int i = 1; i < n-1; ++i){
      apbi = 2.0*(i+1) + apb;
      a[i] = a2b2/((apbi-2.0)*apbi);
      b[i] = sqrt(4.0*(i+1)*(i+1+alpha)*(i+1+beta)*(i+1+apb)/
      ((apbi*apbi-1)*apbi*apbi));
    }

    apbi   = 2.0*n + apb;
    a[n-1] = a2b2/((apbi-2.0)*apbi);

    // find eigenvalues
    TriQL(n, a, &b[0]);

    b.clear() ;
    
    return;
  }


  /** \brief QL algorithm for symmetric tridiagonal matrix

    This subroutine is a translation of an algol procedure,
    num. math. \b 12, 377-383(1968) by martin and wilkinson, as modified
    in num. math. \b 15, 450(1970) by dubrulle.  Handbook for
    auto. comp., vol.ii-linear algebra, 241-248(1971).  This is a
    modified version from numerical recipes.

    This subroutine finds the eigenvalues and first components of the
    eigenvectors of a symmetric tridiagonal matrix by the implicit QL
    method.

    on input:
    - n is the order of the matrix;
    - d contains the diagonal elements of the input matrix;
    - e contains the subdiagonal elements of the input matrix
    in its first n-1 positions. e(n) is arbitrary;

    on output:

    - d contains the eigenvalues in ascending order.
    - e has been destroyed;
  */
  template<typename BaseType>
  void Jacobi<BaseType>::Utility::TriQL(const unsigned int n, Jacobi<BaseType>::UnmanagedArray d, Jacobi<BaseType>::UnmanagedArray e){
    int m,l,iter,i,k;
    BaseType s,r,p,g,f,dd,c,b;
    
    for (l=0;l<n;l++) {
      iter=0;
      do {
        for (m=l;m<n-1;m++) {
          dd=fabs(d[m])+fabs(d[m+1]);
          if (fabs(e[m])+dd == dd) break;
        }
        if (m != l) {
          if (iter++ == MAX_ITERATIONS){
            fprintf(stderr,"triQL: Too many iterations in TQLI");
            exit(1);
          }
          g=(d[l+1]-d[l])/(2.0*e[l]);
          r=sqrt((g*g)+1.0);
          g=d[m]-d[l]+e[l]/(g+sign(r,g));
          s=c=1.0;
          p=0.0;
          for (i=m-1;i>=l;i--) {
            f=s*e[i];
            b=c*e[i];
            if (fabs(f) >= fabs(g)) {
              c=g/f;
              r=sqrt((c*c)+1.0);
              e[i+1]=f*r;
              c *= (s=1.0/r);
            } else {
              s=f/g;
              r=sqrt((s*s)+1.0);
              e[i+1]=g*r;
              s *= (c=1.0/r);
            }
            g=d[i+1]-p;
            r=(d[i]-g)*s+2.0*c*b;
            p=s*r;
            d[i+1]=g+p;
            g=c*r-b;
          }
          d[l]=d[l]-p;
          e[l]=g;
          e[m]=0.0;
        }
      } while (m != l);
    }

    // order eigenvalues
    for(i = 0; i < n-1; ++i){
      k = i;
      p = d[i];
      for(l = i+1; l < n; ++l)
      if (d[l] < p) {
        k = l;
        p = d[l];
      }
      d[k] = d[i];
      d[i] = p;
    }
  }

  // Explicit Instantiation
  template class Jacobi<double>::Utility ;

} // end of namespace





